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We quantize graphs (networks) which consist of a finite number of bonds and vertices. We show 
that the spectral statistics of fully connected graphs is well reproduced by random matrix theory. 
We also define a classical phase space for the graphs, where the dynamics is mixing and the periodic 
orbits proliferate exponentially. An exact trace formula for the quantum spectrum is developed in 
terms of the same periodic orbits, and it is used to investigate the origin of the connection between 
random matrix theory and the underlying chaotic classical dynamics. Being an exact theory, and 
due to its relative simplicity, it offers new insights into this problem which is at the fore-front of the 
research in Quantum Chaos and related fields. 
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■ Quantum graphs (networks) of one-dimensional wires connected at nodes were introduced already more than half a 
century ago to model physical systems. To the best of our knowledge, they appeared for the first time in connection 
with free electron models of organic molecules Jl|-^1 • The molecules were visualized as a set of atoms at fixed locations 
connected by bond paths, along which the electrons obey a one-dimensional Schrodinger equation with an appropriate 
potential. In recent years the interest in quantum graphs has been revived in many areas of physics, and in particular 
in the context of condensed matter physics. Amongst the systems which were successfully modeled by quantum 

""^5 ■ graphs we mention e.g., studies of super-conductivity in granular and artificial materials single-mode acoustic and 
electro-magnetic waveguide networks pO^pdJ , Anderson transition |l2[ and quantum Hall systems with long range 

'-^ ■ potential fL3|j , fracton excitations in fractal structures 14 1, and mesoscopic quantum systems ]lq| . The construction 
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of the wave equations for such networks is a topic in its own right. Ruedenberg and Scherr |4J] (see also |8|]), who were 
apparently among the first to address the problem, based their formulation on the analysis of the limit of wires of 
finite thickness. Quantum graphs can be considered as idealizations of physical networks in the limit where the widths 
O ■ of the wires are much smaller than all the other length scales in the problem. Thus, neglecting the lateral size of the 
, wire, i.e., assuming that the propagating waves remain in a single transverse mode, one replaces the corresponding 
partial differential Schrodinger equation by an ordinary differential operator. This can be justified assuming that 
the inter-mode coupling involves a dynamical tunneling and therefore it diminishes exponentially with the decreasing 
wire thickness. Moreover, when no external field is applied, the motion on the bonds is free, and the problem can be 
further reduced to finite matrices |) 16|. Alexander Jj| was probably the first to discuss networks in external magnetic 
fields. 

Quantum graphs attracted the attention of the mathematics community as well. J. P. Roth |l7j was probably the 
first to derive a trace formula for the spectrum of a Laplacian on graphs. Recently the problem came of age in a series 



of mathematical works by Exner and Seba |1J,|19| , Avron |L6 20 1 and Carlson , whose formulation is based on the 
von-Neumann theory of self-adjoint extensions of formal differential operators (see also J22| and references therein). 

In spite of all this activity, the statistical properties of the spectra displayed by quantum graphs were hardly 
investigated in the past |l^,^3|. Our motivation for studying these spectral properties comes from the theory of 
Quantum Chaos which deals with quantum systems exhibiting chaotic motion in the classical limit. One of the main 
observations of this field is that in the extreme case, when classical motion is strongly chaotic, and in the limit h — > 0, 
the statistical properties of spectra are well described by Random Matrix Theory (RMT). At the same time, the 
spectra of quantized integrable systems display Poissonian statistics. An important goal of quantum chaology, is to 
develop a theory, which relates the quantum spectral statistics to the underlying features of classical dynamics. The 
main tool in this endeavour consists of trace formulae which provide an expression for the spectral density in terms of 
classical periodic manifolds - isolated orbits for chaotic systems, and tori for integrable ones. In most cases, only the 
semiclassical approximation for the trace formulae are known and their application is not only hampered by the 
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intrinsic complexity of the set of periodic orbits, but also by the doubts about the ability of the semiclassical trace 
formulae to provide accurate enough basis for the further developments. This motivated the introduction and the 
study of particular "toy" systems where the required periodic-orbit data can be easily accumulated, while at the same 
time the trace formulae are exact rather than semi-classical p5| . Unfortunately, only very few models combine the 
desirable features of both behaving "typically" and being mathematically simple. It is the main purpose of this paper 
to propose quantum graphs as a very convenient and rich class of systems where the above mentioned requirements 
are met satisfactorily. We shall show that for quantum graphs, one can write an exact trace formula, which is based 
on "periodic orbits" in a way which is analogous to the known trace formulae for chaotic systems. Moreover, we shall 
define the corresponding underlying classical dynamics, and write down the relevant Frobenius Perron operator for 
the "phase space" evolution. This analogy will enable us to study further the connection between spectral fluctuations 
and the classical dynamics. Another clear advantage of quantum graphs is the relative ease by which a large number 
of spectral data can be computed. This enables a rather accurate numerical studies of systems and problems for which 
analytical results are lacking or insufficient. 

This paper extends our previous report f23[ | on the spectral properties of quantum graphs, both in detail and depth. 
In particular, we address issues which have to do with the transition which the spectral statistics undergo when the 
connectivity of the graph is altered. This problem is intimately connected with the semiclassical theory of Anderson 
localization. We also add another parameter to the model, which enables us to interpolate between Dirichlct and 
Neumann boundary conditions. We show that this induces a transition between integrable and chaotic dynamics, and 
we study its effect on the spectral statistics. 

The paper is structured in the following way. In Sec. II, the mathematical model is introduced and the main 
definitions are given. In Sec. Ill we are describing three methods of quantizing graphs. The first method is rather 
standard (see e.g., |l6|] ) and it is the most convenient for numerical computations. The two other methods are related 
to the scattering approach to quantization p6| . One of them forms the basis for the development of the exact trace 
formula, and it singles out a unitary matrix of dimension 2B x 2B where B is the number of bonds on the graph. 
This matrix is the main building block of our theory, and we refer to it as the "bond scattering matrix" Sb- An 
alternative quantization condition is achieved in terms of the "vertex scattering matrix" Sv ■ It is of dimension V x V 
where V is the number of vertices, and it describes the transport through a system in which each of the vertices is 
attached to a conducting wire. In Section IV we present the trace formula for the quantum graph, and also express 
the spectral C function as a sum over composite periodic orbits. Section IV terminates with the introduction and the 
discussion of the underlying classical system. In Section V, the statistical properties of the eigenphase spectrum of 
the bond scattering matrix Sb and of the energy (or wavenumber) spectrum are analyzed and compared with the 
predictions of RMT and of periodic orbits theory. In Section VI, we analyze two families of graphs which are not 
uniformly connected. The resulting spectral statistics deviate from the expectations of RMT, and we explain these 
deviations using periodic orbit theory. Within this study, we investigate also the localization/delocalization transition 
experienced by graphs as a function of the connectivity. Our conclusions are summarized in the last section (Section 



In this section we shall present and discuss the Schrodinger operator for graphs. We start with a few definitions. 
Graphs consist of V vertices connected by B bonds (or edges). The valency Vi of a vertex i is the number of bonds 
meeting at that vertex. The graph is called v-regular if all the vertices have the same valency v. When the vertices 
i and j are connected, we denote the connecting bond by b — {i,j). The same bond can also be referred to as 
b = (Min(i, j), Max(i, j)) or b = (Max(i, j), Min(i, j)) whenever we need to assign a direction to the bond. Several 
bonds connecting the same two vertices are called multiple bonds and the corresponding graph is called a multi-graph. 
Finally, a bond with coinciding endpoints is called a loop. In what follows, unless explicitly specified, we shall consider 
graphs without multiple bonds or loops. Moreover, we shall treat only connected graphs. 

Associated to every graph is its connectivity (adjacency) matrix Ci_j. It is a square matrix of size V whose matrix 
elements Ci.j are given in the following way 



For loop-less graphs the diagonal elements of C are zero. The connectivity matrix of connected graphs cannot 
be written as a block diagonal matrix. The valency of a vertex is given in terms of the connectivity matrix, by 



VII). 



II. QUANTUM GRAPHS: DEFINITIONS 




(1) 
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Vi = J2j=iCij and the total number of bonds is B = ^J2ij=i^iJ- As wm be shown bellow (Sec. IV. C), the 
topological characterization of the graph which was given above, is sufficient for the study of "classical dynamics" on 
graphs. 

For the quantum description we assign to each bond b — a coordinate Xi t j which indicates the position along 
the bond. Xij takes the value at the vertex i and the value Lij = Ljj at the vertex j while Xj^ is zero at j and 
Lij at i. We have thus defined the length matrix Li^ with matrix elements different from zero, whenever Cjj =/= 
and Li t j — Ljj for b = 1, ...,B. The derivations presented in the sequel are valid for any choice of the lengths Li.j. 
However, in some applications we would avoid non generic degeneracies by assuming that the are rationally 
independent. The mean length is defined by (L) = -g Efc=i ^b- 

The wavefunction W is a B— component vector and will be written as (vl/f^ (xb 1 ), ^b 2 ( x b 2 ), ^b B ( x b B )) T where the 
set {bi]f =1 consists of B different bonds. We will call ^b{xb) the component of ^ on the bond b. The bond coordinates 
Xb were defined above. When there is no danger of confusion, we shall use the shorthand notation \&(,(x) for \&f,(x&) 
and it is understood that x is the coordinate on the bond b to which the component ^b refers. 

The Schrodinger operator (with h = 2m = 1) is defined on the graph in the following way PUHl: On each bond 6, 
the component 'J'b of the total wave function \1/ is a solution of the one - dimensional equation 

(-^-^b) = k 2 V b (x), b=(i,j). (2) 

We included a "magnetic vector potential" Af, (with 9fJe(A{,) ^ and — —A^-) which breaks the time reversal 
symmetry. In most applications we shall assume that all the A^s are equal and the bond index will be dropped. 

The wave function must satisfy boundary conditions at the vertices, which ensure continuity (uniqueness) and 
current conservation. The imposition of these boundary conditions guarantees that the resulting Schrodinger operator 
is self-adjoint. The continuity condition requires that at each vertex i, the wave function assume a value denoted by <fi 
which is independent of the bond from where the vertex is approached. Current conservation imposes a condition on 
the derivatives of the wave function at the vertices. By assuming that {h}f =1 — {bi}f =1 the conditions are explicitly 
specified in the following way. For every i = 1, • • • , V: 

• Continuity : 

^i,j(x)\ x -o = Vi, ^i,j( x )\x=Lij = V] f° r a U i < i and C itj ^ 

• Current conservation: 

The parameters Xi are free parameters which determine the boundary conditions. In many applications we shall 
assume that the Xi are all equal, and in such cases the vertex index will be dropped. In the case when Vi = 2, the 
matching conditions can be represented by a (^-function potential of strength Xi . By analogy, we shall refer to the Xi 
as the vertex scattering potential. In the sequel, we shall always assume that > and will consider the domain 
k 2 > (which excludes states bounded at a single vertex). The special case of zero Aj's, corresponds to Neumann 
boundary conditions. Dirichlet boundary conditions are introduced when all the A; = oo. This implies cfi = for all 
i, thus turning the graph into a union of non interacting bonds. A finite value of Xi introduces a new length scale. 
It is natural therefore, to interpret it in physical terms as a representation of a local impurity or an external fields 
JT8| . [l9|j27| . We finally note that the above model can be considered as a generalization of the Kronig-Penney model 
to a multiply connected, yet one dimensional manifold. 



III. THE SPECTRUM OF QUANTUM GRAPHS 

The set of boundary conditions (||), discussed in the previous section, ensures that the Schrodinger operator (||) 
is self-adjoint, and hence the existence of an unbounded, discrete spectrum {fc^}. In the following three subsections 
we shall introduce three different approaches which can be used for the calculation of the wavenumbers spectrum 
{k n }. These approaches complement each other and enable us to address various aspects of the quantum graphs 
using different points of view. 
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A. The Vertex Secular Equation 



The eigenf unctions of the graph are completely determined by their values at the vertices {</>i} i=1 - The quantization 
condition which is to be derived here, specifies the values of k for which a non trivial set of {</?i}ili' s can be found 

The wave function '5 is constructed from B components which correspond to the various bonds. At any bond 
b = the component ^j, can be written in terms of its values on the vertices i and j as 



sin kLij 



(jPi sin [k(L it j - a;)] + ipje %Ai >* Li >J sinfca;) C hj , i < j. 



$ has, by construction, a unique value on the vertices and satisfies the Schrodinger equation 
conservation condition (^) leads to 



(4) 

The current 



sin(feiij) 



V- 

Z / Hi 



kC. 



sin(kLij) 



Vi 



This is a set of linear homogeneous equations for the <p's which has a nontrivial solution when 

det (hi t j(k, A)) = , 



where 



hi j\k, A) 



C iJ e- i ^.^^(sin(fci iJ ))- 1 , z^j . 



The terms /i; )m = ft^ ( , /ly and /i m , m in (Q) diverge when fc is an integer multiple of j^- 
rcctified by replacing the diverging terms by 



hi tm (k,A) = hm,i(k,A) = 
h t j(k,A) = - Ciji cot(kLij>) - 

(k,A) = - C jl>m cot(feXj/ |Wt ) 



(5) 
(6) 
(7) 

This can be easily 
(8) 



k ' 

Am 



The secular equations (0-||) for the quantized graph can be solved numerically to provide an arbitrarily large 
sequence of eigenvalues {k n }. 

As we said previously, the effect of the Dirichlet boundary condition on all the vertices (Aj = oo, Vi) is to disconnect 
the bonds. The eigenfunctions in this case have a simple structure: they vanish on all the bonds except on one of the 
bonds, b, where 



* 6 = -=■ sin(-^— ), JfcW = -f- for all n b > 



(9) 



for all b. The spectrum is the union of the individual spectra, and when the lengths are rationally independent, 
the resulting spectrum displays some Poissonian features. When the \i are large but finite there is always a small 
probability of "leaking" of the wave functions through the vertex scattering potential. As k increases, however, all 
intermediate boundary conditions converge to the Neumann limit. 
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B. Scattering Approach - The bonds Sb Matrix 



The quantization of graphs can be accomplished in a different way which although less efficient from the numerical 
point of view, provides us with a natural and convenient starting point for the construction of the trace formula. It is 
an example of the scattering approach to quantization |p6| . Another variant of this method will be presented in the 
following subsection. 

We first introduce the scattering matrix related to a single vertex. This is done by solving an auxiliary problem of 
a single vertex i, say, with Vi emanating bonds which extend to infinity. The wave- function has components on 
all the bonds bj = (i, j), j = 1, • • • , Vi , which emerge from the i vertex (note that we enumerate and denote the bonds 
in the auxiliary problem by their analogues on the original graph). tpW can be written as a linear combination of 
functions ty^'i) which are solutions for the case where there is an incoming wave entering i from bj = (i, j) and outgoing 
waves from i to all bonds bji (including j = j' which correspond to the reflected part). is a dimensional 

vector with components (xy) for all 1 < f < «j , 

\&f; j) ( Xj ,) =S jtj ,e- ikx i' +iAi >J x i' +af] j ,e ikx i' +iA *-i' x i' . (10) 

Here the Xji are the distances from the vertex i along the bonds (i,f), and crj 8 j, is the v% x m scattering matrix, which 
provides a transformation between the incoming and the outgoing waves at the vertex i. The matching conditions 
(||) at the vertex (xj = 0) together with (|l0|), can be used to determine o,'L: 



af], = (-5 3 . r + + G ""'^ OjCij, , = 2 arctan \ . (11) 

For the Dirichlet boundary conditions we get crj J, = —&j.j< which indicates total reflection. For the Neumann boundary 

condition cr^j, = —5j,j> + ^ which is independent of k. For any intermediate boundary condition, the parameter that 
controls the scattering process is the k dependent parameter 

A - - is ■ < 12 » 

The scattering matrix approaches the Neumann expression as k — * oo. Note that in all the non trivial cases (vi > 2), 
back-scattering (j — j') is singled out both in sign and in magnitude: has always a negative real part, and the 

reflection probability |cj J ]| 2 approaches 1 as the valency Vi increases. One can easily check that is a symmetric 
unitary matrix, ensuring flux conservation and time reversal symmetry at the vertex. For the Neumann boundary 
conditions is a real orthogonal matrix. 

We now write the general expression for an eigenfunction of the quantum graph in terms of its components on the 
bonds b = We write the same bond wave function in two ways. First, we use the standard notation introduced 

before and call it ^>ij(xj) where Xj is the distance from i. The second representation employs the "time reversed" 
notation where the wave function is denoted by ^ j,i{xi) and Xi is the distance from j. The general expressions for 
the wave function in the two representations read 



Vijfa) = a id e i{ - k+Ai >t )x t + b id e i( - k+Ai ^ (13) 
vl> ■ (t ) — n .p i (~ k + A i.*) x i 4- h ■ . P t ( k + A 3,') x i 

= a ■ i e i ^ A s^ x i e i {- k + A i,i) L i,i + b ■ i e i ( k+A J-^ Li -j e i{ -~ k ~ A] ^ X] 



The two representations describe the same function. This gives, 

h,j = a j>i e- ikL *-i- iA ^ L *-i, b jti = a i>j e- ikLi -i- iA ^ Li -i . (14) 

In other words, but for a phase factor, the outgoing wave from the vertex i in the direction j is identical to the 
incoming wave at j coming from i. The incoming and outgoing components of the wavefunction impinging on the 
«'th vertex satisfy 

b h3 = £°S'°ij' ' ( 15 ) 
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( |l4|jl5| ) can be combined to a set of 2B homogeneous linear equations for the coefficients a^j which describe the wave 
function on each of the bonds, a^j is the amplitude for propagation from i to j vertex along the bond b = (i, j), while, 
dj t i is the amplitude for the propagation in the "time reversed" direction (i.e. from j to i) along the same bond. This 
distinction corresponds to assigning directions to the bonds, so that b and b are considered as different entities. We 
thus see that the present approach for quantizing the graph singles out the description in terms of directed bonds as 
the natural setup. 

The condition for a non trivial solution for the 2B dimensional vector (a? , • • • , ar ,a*r , • • • , at ) T gives the secular 

j , v 01 ' ' Ob ° 1 b B ' ° 

equation for the total graph, of the form []2l| 

Cb(A) =det[I-S B (k,A)] = 0. (16) 

Here, the "bond scattering matrix" Ssik, A) = D(k; A)T is a unitary matrix in the 2B dimensional space of directed 
bonds, ft is a product of a diagonal unitary matrix D(k, A) which depends on the metric properties of the graph, and 
a unitary matrix T which depends on the connectivity and on the vertex scattering potentials. 

Dij,i, r (k,A) = ^,^i,j'e ittw+Wyi « ;with L tj = L jt and A l} = -A }i (17) 

T- — A C C rr^ 

The matrix elements of T assign an amplitude to a transition from one directed bond to another. Such a transition 
can occur only if the directed bonds are connected, that is, one is incoming and the other is outgoing from the same 
vertex. The phase and magnitude of the amplitude is given by the corresponding matrix element of the single vertex 
scattering matrix. ^From (|16|) it follows that k n belongs to the wavenumber spectrum if and only if SB(k n , A) has an 
eigenvalue +1. As no approximations were made at any step of the derivation, this quantization condition is exact. 
Furthermore, this gives a constructive method to obtain not only the eigen-energies, but also the wave function, in 
terms of the eigenvector of SB(k n , A) with the eigenvalue +1 |2q] . 

The "bond scattering matrix" Sb cannot be associated with an actual scattering system in the usual sense of 
scattering theory. Nevertheless we shall keep referring to it as a scattering matrix, since it yields a quantization 
condition which is of the standard form in the scattering approach. 

We finally comment that we use the letter £b to denote the secular function (|l6|), because it can be cast in a form 
which is reminiscent of the Riemann-Siegel expression for the Riemann £ function on the critical line. This will be 
shown in chapter IV. 



C. Scattering Approach - The vertex SV Matrix 

The vertex scattering matrix Sv is obtained by converting the graph of interest into a proper scattering system. 
This is done by attaching a lead which is extended to infinity at each of the graph vertices. A scattering solution with 
an incoming wave only in the lead I, and outgoing waves on all the leads can be written in the following way. On the 
external leads, 

^ l) (x) = 5 hl e-* kx + (S v ) h ie lkx . (18) 

On the B internal bonds , 



" 1,3 1 



sin [(Lij - x)k] + ip\ l) e- Ui -^ sinkx) , i < j. (19) 

Sill rvlJi j \ / 



By applying the continuity and current conservation conditions (g) at all the vertices, we get 

<Pi = Si,l + (Sv)i,l 
(S u + {S v )i,i)i + EL h i,j(k)v>T = 



(20) 



where h(k) is the secular matrix defined in (0,||). Combining the above two equations we finally get for Sy 

S v = (U + hik))- 1 ^- h{k)). (21) 

where / is the V x V unit matrix. 



G 



Sy is unitary since h(k) is hermitian, which ensures current conservation. The graph spectrum can be identified as 
the set of wavenumbers for which Sy has 1 as an eigenvalue. This corresponds to a solution where no current flows in 
the leads so that the continuity equations are satisfied on the internal bonds (see Q ) . 1 is in the spectrum of Sv if 

Cv(k) = det [I — Sv] = < — ► 2 V det [il + h{k)]~ l det h(k) = (22) 

which is satisfied once det h(k) = 0. This is identical with the condition (JgJ) which was derived in subsection III. A. 

Variations on the same theme can be obtained by considering graph scattering systems where leads are attached 
to an arbitrary set of L vertices {ty}/ =1 , with 1 < L < V. The L x L scattering matrix Sv has to be modified in the 
following way 

S v = (h(k) + iW T W)~ 1 W T - I (23) 

where = is the LxV leads - vertices coupling matrix. In the case that we examined previously with L = V, 
W = I. 

The matrices Sv will not be studied any further in this work, and the derivation above was given for the sake 
of completeness. The Sv scattering matrix corresponds to proper scattering problems, and can be used to model 
experimental systems, such as e.g., conductance of mesoscopic microdots. We shall study the statistics of conductance 
fluctuations, based on the vertex scattering matrices Sv in a separate publication |2S|| . 

IV. PERIODIC ORBITS, THE TRACE FORMULA AND CLASSICAL DYNAMICS ON GRAPHS 

In this section we derive an expression for the quantal density of states in terms of periodic orbits on the graph. A 
trace formula for the Laplacian on graphs was first presented by J. P. Roth Jl7|] - Our result (see also pi ) generalizes 
Roth's expression in several ways, and it is derived by means of a different approach. The key element in this theory 
is the concept of a periodic orbit on the graph, which we shall introduce at this point. 

An orbit on the graph is an itinerary (finite or infinite) of successively connected directed bonds (ii, 12), («2> ^3)) " • '■ 
For graphs without loops or multiple bonds, this is uniquely defined by the sequence of vertices *i, ia, - * - with i m S 
[1, V] and Ci m> i m+1 = 1 for all m. An orbit is periodic with period n if for all k, (i n +ki hi+k+i) — (ik, ik+i)- The code 
of a periodic orbit of period n is the sequence of n vertices i\, ■ ■ ■ ,i n and the orbit consists of the bonds (i m , i m +i) 
(with the identification i m +n = im)- I n this way, any cyclic permutation of the code defines the same periodic orbit. 

The periodic orbits (PO) can be classified in the following way: 

• Irreducible periodic orbits - PO's which do not intersect themselves so that any vertex label in the code can 
appear at most once. Since the graphs are finite, the maximum period of irreducible PO's is V. To each 
irreducible PO corresponds its time reversed partner whose code is read in the reverse order. The only code 
which is both irreducible and conjugate to itself under time reversal is the code corresponding to PO's of period 
2. 

• Reducible periodic orbits - PO's whose code is constructed by inserting the code of any number of irreducible 
PO's at any position which is consistent with the connectivity matrix. All the PO's of period n > V are 
reducible. 

• Primitive periodic orbits - PO's whose code cannot be written down as a repetition of a shorter code. 

After these preliminaries, we are set to derive the trace formula for the graphs. Once this is done, we shall show 
that one can define classical dynamics on the graph, and that the periodic orbits on the graph are indeed the analogue 
of the periodic orbits of hyperbolic classical Hamiltonian systems. 

A. The Trace Formula 

The starting point for the derivation is the secular equation (|l6|). The function Cb(&) is a complex valued function. 
It will be convenient to write it as a real amplitude times a phase factor. Denoting the eigenvalues of Ssik, A) by 
e i0,O) for I = 1,...,2B we get 
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C s (fc)=exp(^e(fc))2 2B nsin^ (24) 



where 0(fc) is 



2D 



Q(k) = ~log(det(-S(fc, A)) = V 9rfk, A) - 2Bn 
i 

1=1 

v 

= kC + (B - V>- 2^arctan(A l ) . (25) 



Here £ = 2 z)2b=i Lb is twice the total length of the graph, and the Aj depend on the boundary conditions as defined 
in (^|). Notice that the parameters Ab do not appear in the above expression. This is because the contributions of 
time reversed bonds are canceled pairwise. 

The last product in ( pi] ) is real on the real k axis. Therefore, the imaginary part of its logarithmic derivative is a 
sum of delta distributions located where Cs(fc) = 0. Using the expansion 

logdet(7 - S B (k)) = - V -trSg(fc) (26) 

' n 

n— 1 

we obtain the following expression for the density of states 

d ^ = ^^T e ( fc ) + -lim^o^m-^- £ ~tiS%(k + ie) . (27) 
2ir ak 7r ak — ' n 

n— 1 

The first term on the right hand side of ( p7[ ) corresponds to the smooth spectral density while the second one provides 
the fluctuating part. 

The spectral counting function N(k) is given by 

N(k) = [ d(k')dk'. (28) 
Jo 

From (|27j) we have : 

oo 

N(k) = N(k) + -9fm V -tr(S B (k)) n (29) 

TV < ^ 71 



n— 1 



where 



N(k) - N(0) = A- [0(A) - 9(0)] = !f + \ - - arctan(A l ) . (30) 

Z7T Z7T Z 7T ' — ' 

i=l 

This is the smooth part of the spectral counting function. The leading term involves the "volume" C of our system, 
and it is independent on the boundary condition A.;. The next two terms are due to the scattering potentials on the 
vertices. The contribution is minus the value of the third term at k = 0. For the Neumann boundary condition, the 
limit A; — > should be taken after the value k = is substituted. For large k the last term is inversely proportional 
to k if Xi 7^ 0. Hence, the mean level density d = dkN(k) is essentially constant, reflecting the fact that the graph 
is one dimensional. For the Neumann boundary conditions (A, = 0), d is independent of the wave number k and 
N(0) = 1/2 and 

N(k) = f + l- (31) 
Ztt Z 

The oscillatory part of the counting function is expressed in terms of tr(SB(k)) n . Using the definitions jT7| ) and 
S = DT one can obtain the tr(5e(fc))™ directly as sums over n— periodic orbits on the graph: 

tr(S B (k)) n = n p A r p e l{k ^ +<bp)r e l ^ +p ^ r (32) 
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where the sum is over the set V n of primitive PO's whose period n p is a divisor of n, with r = n/n p . l p — J^bep ^ b 
is the length of the periodic orbit. <5> p = J^bep ^bAb is the "magnetic flux" through the orbit. If all the parameters 
Ab have the same absolute size A we can write $ p = Ab p , where b p is the directed length of the orbit. fi p is the 
number of vertices (with m > 2) where back scattering occurs. At the other v p vertices on the PO the scattering is 
not backwards. The number of back scatters from vertices with Vi — 1 is n p — (p p + v p ) > 0. The amplitudes A p are 
given by 



(i- 



n 



v t (l + iA t ) 



{l + iA s ) 

where 7 p plays the role of the Lyapunov exponent. When A.; =/= 0, the phase of each term is modified by adding 

2A„ 



(33) 



Hp 

p p (k) = y arctan 



i=l 



Substituting (t 



> « i (l+A|)-2 
in (|27D one gets an exact trace formula 



^ 2 

i=l 



n p -(/Lip+Kp) 

E < 

i=l 



d(*0 = 75- + - 



1^ 



A, 



2tt 



7T 



1 lp cos { r {kl p 

IT 



(34) 



(35) 



p p ir + p p (k))) 



e 2 



The above formula bears a striking formal similarity to the well known exact Selberg trace formula |2£] for modular 
domains on Riemann surfaces with constant negative curvature, to the Riemann Weyl relation for the Ricmann zeros 
on the critical line and to the semi-classical Gutzwiller trace formula for chaotic Hamiltonian systems |24|. As we 
shall show in the sequel, the analogy between PO's on the graph and periodic orbits in dynamical systems follows 
naturally from the classical dynamics which we associate with the graph. This analogy is strengthened by further 
evidence: The number of n— PO's on the graph is itrC n , where C is the connectivity matrix. Since its largest 
eigenvalue T c is bounded between the minimum and the maximum valency i.e. mhiVi < T c < maxw^, periodic orbits 
proliferate exponentially with topological entropy ps logr c . The amplitudes A p which play the role of the stability 
amplitudes, decrease exponentially with n but not enough to make the series for d(k) absolutely convergent (positive 
entropy barrier). Finally, /i p is the analogue of the Maslov index. Its origin is topological, and it counts the number 
of non trivial back scatters along the PO. This can be expressed as the number of sequences of strings of the type 
• • • i a , if,, i a ■ ■ ■ (with vi, > 2) which appear in the code of the PO. 

The distinguishing feature of the graph trace formula is the structure of the spectrum of lengths l p of the periodic 
orbits which appear in (|35|). In contrast with the other systems mentioned above, the lengths are constructed as 
linear combinations (with integer coefficients) of the elementary bond lengths Lf,. Hence, the lengths spectrum is 
characterized by a high degree of degeneracy. The degeneracy is not too important for PO's of period n < V. As n 
increases, it becomes progressively dominant, and it is the main feature of the length spectrum for periodic orbits with 
n > 2B. The effects of this degeneracy are most apparent in the study of spectral statistics which will be discussed 
in the next chapter. 



B. Periodic orbits expression for the spectral £ function 

The spectral density ( |35| ) is not convenient to deal with, since it is not a proper function but, rather, a distribution. 
In this section we would like to cast the information which is stored in ( |35| ) in a different form, and express it using 
periodic orbits. 

The first method for achieving this goal is based on the identity (^6|), where the periodic orbit expression ([32j) for 
trSg is substituted, and the summation over the repetitions is carried out explicitly. One gets 

C fl (fc) =JJ(l-tp) ; t p =e-^e i{kl ^ +Ah ^S^ +p ^ k)) . (36) 
p 

This expression as a product over primitive PO's justifies the letter £ by which the secular function is denoted. This 
product does not converge on the real k axis because the number of primitive PO's proliferates exponentially, with a 
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topological entropy which approximately equals the mean Lyapunov exponent 7. Hence, the product (|36|) converges 
in the absolute sense only for Smfc > j^y- This is the "entropy barrier" for the ( B function. As a matter of fact, 
the formal manipulations which were used to derive ( |36| ) are strictly justified beyond the entropy barrier. 

A converging, well behaved expression for the (b function on the real k axis can be obtained, and it is the analogue 
of the Ricmann Siegel expression for the Riemann £ on the critical line. To this end, we first note that the Cb(^) 
function (|l6|), is the characteristic polynomial of the matrix S B {k) 



2B 



<; B (z,k) = det(zl - S B {k)) =Y, ai (k)z l (37) 

;=o 

evaluated at z = 1. The coefficients a n (k) satisfy an inversive symmetry relation which follows from the unitarity of 
the scattering matrix S B (k) 

a 2B _ l (k) = e ie ^a*(k). (38) 
Utilizing relation (|3^ ) one may rewrite < ^T\) in a more convenient form 

ynJi;,, =2e i ^Rele- i ^r 1 (y ai + l llLS I }. CW) 




We would like to emphasize that this form of the secular function ( B (k) is due to the unitarity of S B (k). Thus, the 
removal of the contributions of terms with I > B in (|3^) is not only a practical saving of numerical effort, but also an 
expression of a basic property of the system. 

Let us now consider the S B (k) matrix at a given wavenumber k. Its spectrum is on the unit circle i.e z = exp(i9) 
and ( B (e l9 , k) becomes a function of 9 which depends parametrically on k. Yet, the secular function ( ^9|) is not real 
for real 9. It is useful to define another function which is real on the real 9 line, and vanishes at 9 = 9\. Thus, we 
introduce the spectral determinant which is obtained from (|37]) by extracting a phase factor 

Z(9,k) = e l 2 Cs(e,&) = 2 2B [[stn [ 

1=1 ^ 

= 2Re je-<^^ (|>- + \a B e^ j . (40) 

The function Z(9, k) serves a dual purpose. Setting 9 = it is a real secular equation for the graph. For a fixed 
value of k it is the secular function for the spectrum of S B (k). 

To write Z{9, k) in terms of periodic orbits it is useful to recall Newton's identities (2(| 

at = -y ^trS' + ^a„trS'-™^ , l = l,...,2B. (41) 



The explicit dependence of the a/ on trS n takes the form |2q| 

ai=-jtrS l B (42) 

— 7 n™ " ; ^ T ^B ltr ^B 1 toS B 2 ' ' ' trS*^ 1 1 n tiSg 
1 — » lli=l l i 

where the summation is over all vectors ^ with integer entries such that I > h > 1% > • • • > In > 1- By substituting 
further from fl3^ ) the trS 1 ^ we get an expression of the a; , for I < B in terms of composite orbits . 

The form obtained above for the spectral £ function is reminiscent of the Riemann-Siegel approximation for the 
Riemann £ on the critical line. Here, however, it is an exact expression, and because of the fact that the density of 
states is constant, the number of terms appearing in the sum is independent of k. 
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C. Classical evolution 



We introduced above the concept of orbits on the graph as strings of vertex labels whose ordering obeyed the 
required connectivity. This is a finite coding which is governed by a Markovian grammar provided by the connectivity 
matrix. In this sense, the symbolic dynamics on the graph is Bernoulli. So far, the orbits were used and discussed 
as formal symbols devoid of a dynamical origin. In the present section we shall introduce the classical dynamics 
which can be associated with a graph, and which complements the quantum dynamics on the graph. The classical 
dynamics makes use of the representation of the graph in terms of directed bonds. In the following sections we shall 
use the label b to refer to the directed bonds, so that b = 1, • ■ • , 2B and we shall denote the time reversed pairs by 
b = ; b = 

We consider a classical particle which moves freely as long as it is on a bond. The vertices are singular points, 
and it is not possible to write down the analogue of the Newton equations at the vertices. Instead, one can employ 
a Liouvillian approach, based on the study of the evolution of phase-space densities. The phase space evolution 
operator assigns transition probabilities between phase space points, for which a quantum analogue can be found. 
The phase-space description will be constructed on a Poincare section which is defined in the following way. Crossing 
of the section is registered as the particle encounters a vertex, thus the "coordinate" on the section is the vertex label. 
The corresponding "momentum" is the direction in which the particle moves when it emerges from the vertex. This 
is completely specified by the label of the next vertex to be encountered. In other words, 

position 1 J vertex index 

momentum J 1 next index 

The set of all possible vertices and directions is equivalent to the set of 2B directed bonds. The evolution on this 

(i) 

Poincare section is well defined once we postulate the transition probabilities P^L,^ between the directed bonds 
b = and b' — (i,f)- To make the connection with the quantum description, we adopt the quantum transition 
probabilities, expressed as the absolute squares of the Sb matrix elements 



p(0 



(l + e -**) 
-hi' + 



(43) 



The phases oji are given in ([11]). For the two extreme cases, corresponding to Neumann and Dirichlet boundary 
conditions, (E3|) results 



Pjtif = (-hf + |) > Neumann (44) 
= 5j ji , Dirichlet. 

(i) 

The transition probability Pj^^r for the Dirichlet case, admits the following physical interpretation. The particle 
is confined to the bond where it started and thus the phase space is divided to non-overlapping ergodic components 
(«"tori"). For all other boundary conditions, the graph is dynamically connected. 

The classical evolution (Frobenius Perron) operator Ub t b' between the bonds b, b' reads 

Uij.nm = dj^Pijl+jm (45) 

U does not involve any metric information on the graph, and for Dirichlet or Neumann boundary conditions U is 
independent of k. 

If pb(t) denotes the probability to occupy the bond b at the (topological) time t, then we can write down a Markovian 
Master equation for the classical density: 

p b (t + l) = J2 U b,b'Pb'(t)- (46) 

b' 

The unitarity of the graph scattering matrix Sb guarantees J^6=i Ub,b' — 1 an d < Ub.b' < lj so that the probability 
that the particle is on any bond is conserved during the evolution. The spectrum of U is restricted to the interior 
of the unit circle and v\ = 1 is always an eigenvalue with the corresponding eigenvector |1) = (1, 1, .... 1) T . In 
most cases, the eigenvalue 1 is the only eigenvalue on the unit circle. Then, the evolution is ergodic since any initial 
density will evolve to the eigenvector |1) which corresponds to a uniform distribution (equilibrium). The rate at which 
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equilibrium is approached is determined by the gap between the next largest eigenvalue and 1. However, there are 
some non generic cases, such as e.g., bipartite graphs when —1 belongs to the spectrum. In this case the asymptotic 
distribution is not stationary (see for example Section VII). If 1 is the only eigenvalue on the unit circle we have 

p(t) t=^ |1) (47) 

with a mixing rate In | z^2 1 determined by z/2, the second largest eigenvalue of U. This is characteristic of a classically 
mixing system. 

Of prime importance in the discussion of the relation between the classical and the quantum dynamics are the 
traces u n = tr(f/ n ) which are interpreted as the mean classical probability to perform n- periodic motion. When only 
one eigenvalue v is on the unit circle, one has that u n n^o& 1. Then, we can obtain a classical sum- rule by substituting 
the periodic orbit expansion of u„, 



n p (\A p \ 2 ) n^5tl. (48) 



Each periodic orbit is endowed with a weight |-A P | 2 defined in terms of the stability amplitudes (|33[). It is the 
probability to remain on the orbit. These weights are the counterparts of the stability weights | det(I — M p )| _1 for 
hyperbolic periodic orbits in Hamiltonian systems, where M p is the monodromy matrix. Graphs, however, are one 
dimensional and the motion on the bonds is simple and stable. Ergodic (mixing) dynamics is generated because 
at each vertex a (Markovian) choice of one out of v directions is made. Thus, chaos on graphs originates from the 
multiple connectivity of the ( otherwise linear) system. 

Using the expression ( |4§| ) for u n one can easily write down the complete thermodynamic formalism for the graph. 
Here, we shall only quote the periodic orbit expression for the Ruelle £ function 

Cr{z) = (det(7 - zU))^ 1 = exp [-tr (ln(J - zU))} (49) 

1 



exp 



V — U n 
^ — ' n 



IV 



(i (1 - z n p exp(-n p 7 p )) 



where the product extends over all primitive periodic orbits and we have used the definitions of (p 



V. SPECTRAL STATISTICS 



So far we developed the spectral theory of graphs, pointing out the similarity between quantum graphs and more 
complex quantum systems which display chaotic dynamics in the classical limit. In the present chapter we shall report 
about analytical and numerical results which show that the spectral statistics of these simple systems also follow the 
pattern of more general Hamiltonian systems. Namely, when the classical graph dynamics is mixing, and in the limit 
when the (topological) time needed to reach equilibrium is much shorter than the number of directed bonds, the 
spectral statistics for quantum graphs are very well reproduced by the predictions of Random Matrix Theory. In 
the integrable limit (Dirichlet boundary conditions) the graph spectral statistics is Poisson as is the case in generic 
integrable Hamiltonian systems. The investigation of the universality of spectral fluctuations and deviations thereof is 
especially convenient for graphs because of the transparent and simple spectral theory in terms of PO's, and because 
of the relative numerical ease by which large spectral data bases can be constructed. The parameters which appear in 
the theory can be used to study various characteristic spectral transitions: Changing the vertex potential parameters 
Xi (which, for simplicity, we will choose to be the same for all the vertices i.e. A = Xi) induces the transition between 
classical integrability to chaos. This is accompanied by a spectral transition from Poisson to RMT like statistics. The 
parameter (again for the sake of clarity we will choose \Ab\ — A for all the bonds) allows to break time reversal 
symmetry, and the mean valency can be used to study the dependence of spectral statistics on connectivity. All these 
questions will be dealt with in the present chapter. 

The spectral theory presented above relied heavily on the bond scattering matrix Sb which provided the foundation 
for the periodic orbits theory and the connection to the classical evolution. The spectral statistics of the eigenphases of 
the Sb matrix are also intimately connected with the spectral statistics of the graph wavenumber spectrum. Therefore, 
we will start our presentation by discussing the Sb matrix spectral statistics. 
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A. The spectral statistics of the Sb matrix 



We consider the SB(k,A) matrices defined in (16, [l?]). Their spectrum consist of 2B points confined to the unit 
circle (eigenphases) . Unitary matrices of the type considered here are commonly studied since they are the quantum 
analogues of classical, area preserving maps. Their spectral fluctuations depend on the nature of the underlying 
classical dynamics |32] l . The quantum analogue of classically integrable maps display Poissonian statistics while in 
the opposite case of classically chaotic maps, the eigenphases statistics conform quite accurately with the results of 
RMT for Dyson's circular ensembles. The ensemble of unitary matrices which will be used for the statistical study 
will be the set of matrices <Ss(fc, A) with fc in the range \k — ko\ < A^/2, where the mid point fco and the interval size 
A/c are to be determined. Since the dimension of the Sb matrices is independent of fc, the mean value fco is important 
only when the boundary conditions are neither Neumann nor Dirichlet. For the intermediate boundary conditions, 



fco sets the mean value of the parameter A (12), and A*, must be chosen such that A does not change appreciably in 
the interval. However, A& must be much larger than the correlation length between the matrices Ss(k, A). One can 
estimate the correlation length by studying the auto-correlation function 



C(X) = -r- / ^tr (4(fc' + f )S B (fc' f )) dk'. (50) 

A fc Jk„-A k /2 2B V 2 2 / 



For the two extreme cases of Neumann and Dirichlet boundary conditions the auto-correlation function (50) can 



be calculated exactly. By writing the scattering matrix as Ss{k) = exp{ikL)Ss{Q) and substituting in (|5C|) we find: 

= \ E ^ = / e lLx P(L)dL ee P( X ) (51) 



i=i 



where P(x) is the Fourier transform of the probability distribution P(L) of the lengths of the bonds. Thus, the 
correlation is inversely proportional to the variance of the distribution of the lengths if,. From now on, we shall 
assume A& 3> var(Ls), which justifies the fc averaging procedure. The ensemble average will be denoted by 



1 



fco+Afc/2 

• dk . (52) 



&k Jfc _A fc /2 

Another way to generate an ensemble of graphs, is by randomizing the length matrix L which contains the lengths 
of the bonds while the connectivity (topology of the graph) is kept constant. This is the disorder approach, which 
will also be applied when called for. 

In the following subsections VI. A. 1-2, we investigate some statistical measures of the eigenphases {8i(k)} of the 
scattering matrix Sb and compare them with the predictions of RMT, and with the results of the periodic 

orbits theory of spectral fluctuations which was originally developed for quantized maps. The two statistics which we 
shall investigate are the spectral form factor and the autocorrelation of the spectral £ function. Explicit expressions 
for these quantities are given by Random Matrix Theory p5| and a semiclassical theory is also available ]3^-]3^] . 



1. The Form Factor 



The Sb matrix for a fixed eigenvalue fc is a unitary matrix with eigenvalues e l6l< - k \ The spectral density of the 
eigenphases then reads 

2B oo 

d(6; fc) ee y 6(6 - 9 t (fc)) = — + — V e~ l6n tr S B (fc)" + c.c (53) 

— ' 27T 27T * — ' 

1=1 n=l 

where the first term on the r.h.s is the smooth density d = while the others describe the fluctuating part. 

The two-point correlations are expressed in terms of the excess probability density R 2 (r) of finding two phases at 
a distance r, where r is measured in units of the mean spacing |j 



-R2<V; fc ) = 
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jtrS'gl 2 >k is the Fourier transform of i?2(r, fco). F° r a Poisson spectrum, 



only 133] , and explicit expressions 



The form factor K(n, 2B) = ^ < 
K(n, 2B) — 1 for all n. RMT predicts that K(n, 2B), depends on the scaled time ys 
for the orthogonal and the unitary circular ensembles are known pEj . 

We computed 55 < (tr^sl 2 >fc f° r we U connected graphs, with various vertex potential parameters. In Fig. 1 we 
show typical examples, calculated for a fully connected pentagon. The results for Neumann boundary conditions show 
quite a good agreement with the predictions of RMT for the Circular ensembles. We shall discuss and explain these 
results in the following paragraphs. 

To begin, consider the data for Neumann boundary conditions and A = or A 7^ (see Fig. 1). The predictions of 
RMT are also shown, and they reproduce quite well the smooth trend of the data in the two cases. The deviations 
from the smooth curves are not statistical, and cannot be ironed out by further averaging. Rather, they are due to 
the fact that the graph is a dynamical system which cannot be described by RMT in all detail. To study this point in 
depth we shall express the form factor in terms of the PO expression (p2[). Assuming Neumann boundary conditions 
for the time being, 



2B V sv 11 Ik 2B 



n ^j{r e i(kl p +Ab p +np, p )r 



(55) 



— ^ n p n pl A r p A r pl eyiY>{iA(rb p -r'b pl )+iTi(r^ p ~r'^ pl )} 



p,p'£V„ 



The k averaging is carried out on such a large interval that the double sum above is restricted to pairs of periodic orbits 
which have exactly the same length. The fact that we choose the lengths of the bonds to be rationally independent 
will enter the considerations which follow in a crucial way. Consider first the domain 1 < n <C 2B. The PO's are 
mostly of the irreducible type, and the length restriction limits the sum to pairs of orbits which are conjugated under 
time reversal. Neglecting the contributions from repetitions and from self tracing orbits we get 

^L(|trSS(fc)| 2 ) fc «^ £ n 2 A 2 4cos *Ab p = ^u n (cos 2 Ab p ) n . (56) 
pev n 



The classical return probability u n approaches 1 as n increases (see (^8|)). However, deviations from unity reflect 
the fact that the classical dynamics reaches the ergodic state only after some time. The deviation which is simplest 
to understand occurs at n = 1. Since there are no classical fixed points (no self connected vertices) on the graph, 
u\ = 0. However in the limit B — > 00, the short time deviations converge to the origin when the scaled form factor 
is studied. Neglecting the short time deviations, we can replace u n by 1, and we see that the remaining expression 
to be evaluated is the classical expectation of cos 2 Ab p over PO's of length n. For A — this factor is identically 1 
and one obtains the leading term of the COE expression for n <C 2B. If A is sufficiently large (cos 2 Ab p ) sa 1/2, and 
one obtains the short time limit of the CUE expression. The transition between the two extreme situations is well 
described by 

(cos 2 ^ p > n «i(e-^^)f +l) . (57) 

This formula is derived by assuming that the total directed length b p of a periodic orbit is a sum of elementary lengths 
with random signs. 

One cannot use the arguments presented above for the range n > B. As n approaches B the degeneracy of the 
length spectrum prevails and for n > 2B all the orbits are degenerate. In other words, the restriction rl p — r'l p i in 
([55|) does not pick up a unique orbit and its time reversed partner, but rather, a group of isometric but distinct orbits. 
Therefore, the interferences of the contributions from the group of all the orbits must be calculated. The relative 
sign of the terms is determined by the "Maslov" index. This can be seen better, once rewriting (p^) in the form (we 
assume for simplicity A = 0) 



J2 n p A r p e™^ r 



(58) 
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where the second sum is over the set V q of PO of the type l p — IbLb with <Zfc — n p . It is clear that the indices of 
different orbits in a family of isometric PO's are correlated. Otherwise, if the Maslov indices are random, one would 
regain the diagonal approximation (^6|) for arbitrarily long times. The correlation between the Maslov indices within 
the family of isometric PO's are the analogue of action correlations in the semiclassical theory of spectral statistics 
p9[ . The dynamical origin of these correlations is not known also for graphs, and it is one of the important open 
problems that should be addressed. 

Graphs with Dirichlct boundary conditions are integrable in the sense explained above. One expects therefore that 
the spectral statistics in this case is Poissonian, which implies 

^(|tr^(fc)| 2 ) fe = l for all n>l. (59) 

In the Dirichlet case, the Sb matrix reduces to a block diagonal form where each bond and its time reversed partner 
are coupled by a 2 x 2 matrix of the form 

... / i(k+A)L b \ 

S {b Hk,A)=^ ei(k _ A)Lb e Q j. (60) 

The spectrum of each block is the pair ±e lkLb , with the corresponding symmetric and antisymmetric eigenvectors 
■^-(1,±1). As a result, we get 

^(|tr^(fc)| 2 ) fe = l + (-l) n for all n > 1 . (61) 

This deviation from the expected Poissonian result is due to the fact that the extra symmetry reduces the Sb matrix 
further into the symmetric and antisymmetric subspaces. The spectrum in each is Poissonian, but when combined 
together, the fact that the eigenvalues in the two spectra differ only by a sign is the reason for the anomaly (|6l|). 
In the sequel we shall remove this feature by considering the smooth form factor obtained by taking the mean of 
successive n values. 

The transition between the two extreme boundary conditions can be affected by using the interpolating boundary 
conditions where A ^ {0,oo}. The relevant parameters are the A defined in (|l2|), and it is expected that the spectral 
statistics make the transition from RMT like to Poisson as these parameters spans their range of values. The exact 
symmetry which prevails in the Dirichlet case, is broken for intermediate values of A. However, the tendency towards 
trapping in a single bond is a dynamical feature, which persists for finite A values, and therefore the even-odd 
staggering of the form factor can be observed also for the intermediate range of A values (see inset in Fig. 1). Since 
the dynamical reason for this effect is clearly understood, we show in Fig. 1 the pairwise averaged form factor, which 
displays clearly the transition from the Poissonian to the RMT (COE and CUE) limit. 

As was mentioned above, the short times deviations of the data from the RMT expectations (see Fig. 1) are real, 
and are due to the deterministic nature of the dynamics induced by the Sb matrices. It is easy to show this explicitly 
for (|tr5 2 | 2 ), since here all the contributions are due to period-2 PO's which are self tracing, and each has its distinct 
length. Using (|32] ) we get 

1 'KWI'), = 2 f<^M\ (62 , 



2B V ay n Ik \ 1 + A 2 

independently of the value of A. This is different from the value -i expected for the CUE and « % expected for the 
COE. 



2. Spectral Z function correlations 

We conclude this section with an analysis of the spectral Z function given by (^) . The statistical properties of this 
function, can be expressed in terms of the statistics of either the eigenphases 6i or the coefficients a; (see (|o|)). Since 
the two sets of variables are functionally related, they are statistically equivalent. In practice, however, one cannot 
check the full spectral distribution, and therefore it is advantageous to study statistical measures which are based on 
other accessible quantities. The measure which was proposed in [p7| and p0| was the autocorrelation function 
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r^M (z(6 + 2 ;k )z*(e-l-,k)) k ^ (63) 



2B 

\ 2 ) P Ml-B) 
k 



1=0 

This statistical measure, depends on higher order correlations of the eigen-phases. Hence, C(rf) tests aspects of the 
eigenvalues distribution which are not accessible by the study of the two-point form factor discussed previously. 

The ensemble averages (\ai\ 2 ^ for Circular Random Matrices were calculated in |3jJ. We shall quote here the 

results for the COE and CUE ensembles: 

= 1, 13 = 2 



An approximate expression for the ^|a/| 2 y was obtained by assuming that the txS n are independent random 

Gaussian variables for n < B p7[ . This approximation is an extension of the diagonal approximation mentioned 
above, and it leads to the following recursion relation 



^i 2 )-}E(i a '-i 2 ) vltl ! B ' ' «*'' 

11=1 

For the calculation of the autocorrelation function (|63|), it is sufficient to obtain the (|a;| 2 ) for I < B. The rest are 
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provided by the inverse symmetry fl38|) which moreover implies that the Fourier components of C{rf) are symmetric 
about B and thus the autocorrelation function is real. 

Using the approximate result (|56[ ) for the form factor ^|trS'g| 2 ^) ~ gnu n we have fjofl 

i 

ai = (\ a i-k\ 2 ) u k, (66) 

k=l 

which should be solved with the initial condition ^|ao| 2 ^ = 1. In ( |66[ ) g = 2(1) for systems with (without) time reversal 

symmetry, and Uk = tiU k is the classical return probability. For systems which display strong mixing, u n = 1, and the 
approximate recursion relations reproduce the RMT result for systems which violate time reversal symmetry (/3 = 2). 
For systems which are invariant under time reversal, one reproduces only the leading term in 4^ of the RMT result 
(/3 = 1). 

We computed numerically the (|a/| 2 ) fe for the completely connected pentagon, subject to Neumann boundary 
conditions, where time reversal symmetry was either respected (A = 0) (Fig. 2a), or violated (A ^ 0) (Fig. 2b). The 
results are displayed in Fig. 2, and they deviate substantially from the RMT predictions (|64|). Note that the values 

of (^trSgl 2 ^ for the same system showed a rather good agreement with RMT (see Fig. 1). The reason for the large 

deviation is clear. No physical system can reproduce the strong mixing condition u n = 1 for all n. Indeed, this is 

the reason why |tr5"|| | 2 ^> deviate from the RMT results for short times. Because of the iterative procedure (|6 

the short time non-generic effects reveal themselves in the higher order coefficients a/ , and this is why this statistical 
measure is much more sensitive to the non universal features of the classical dynamics. The approximate theory 
presented above, includes the correct short time behavior of the system, and therefore it reproduces the main features 
of the numerical data much better than RMT. A quantitative measure for the expected deviation from the RMT 
prediction can be given by the magnitude of the next to the leading eigenvalue of the classical evolution operator U . 
For the system we studied, it is -.25, which is still far from the value expected in the strong mixing limit. 

Let us finally comment on the transition from Poisson to RMT due to variation of the parameters Aj. As it 
was mentioned already in the introduction, for integrable systems, we expect that the spectrum is uncorrelated and 
described by the Poisson ensemble which gives for the (|a/| ) the expression 



H 2 } = ( 7 ] ■ (67) 



2B 
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On the other hand, the diagonal approximation predicts that (|trS'™| 2 ) ~ 2B (see previous section) and thus it 
provides us with the following recursion relation for the coefficients (|a;| 2 ) 

/I l 2 \ 2B (l a /-n| 2 > (czq\ 

n=l 

In Fig. 2a, b we present our numerical results for a fully connected pentagon and for various values of the parameter 
A. Again, we see that our system undergoes a transition from GOE/GUE to Poisson statistics when A increases. 



B. Level Statistics 

The statistical properties of the energy levels (or the wave numbers) spectrum can be derived from the statistics of 
the eigenphases of Sb because of the following reasoning. The wavenumber spectral density can be written as 



2B 



= 1 



dk 



(69) 



One can easily show that 



L < M < L (70) 
dk 



where L m i n ,max denote the minimal or the maximal bond length of the graph. Consider an interval dk about fco so that 
\ d dfc fc ^ / ^ = ^ ^ ^ 7T '■ Since the mean wavenumber spectral density is ^ B the interval dk can accommodate a 
large number of levels when B is large. The wavenumber density in the dk vicinity of fco is 



d(fc;fc ) = i>. (eM + ik-ko)^^ ^jp«<£>£ <MW + (* -k ) (L)) (71) 

For a given k value, the expression on the rhs is the eigenphase density of the unitary matrix 5b (&o), (the I independent 
shift of the phases does not change the distribution of intervals on the circle). Hence, one can read the short range 
statistical properties of the k spectrum, from the results on the statistics of the eigenphases which was discussed in the 
previous section. In the sequel we shall supply numerical data and additional arguments to show that this is indeed 
the case. We shall also compute various statistical measures which are commonly used in the statistical analysis of 
spectral fluctuations of quantum systems whose classical analogue is chaotic. We shall show that the spectrum of the 
quantized graph behaves as a typical member of this set of "quantum-chaotic" systems. 

For the numerical calculation of the spectrum we had used the method described in section III. A. That is, we 
identified the spectrum as the zeros of det h(k, A). The completeness of the spectrum was checked by comparing the 
counting function N(k) with Weyl's law (|30|). An efficient detector of missing or spurious levels is provided by the 
function S n 

5 n = N(k)-N(k) . (72) 

This quantity is expected to fluctuate around zero and a redundant or a missing eigenvalue is accompanied by an 
offset by ±1. Fig. 3 shows a typical plot of S n as a function of n for a graph containing 5 vertices (see inset in Fig. 3) 
and Neumann boundary conditions. S n fluctuates about as expected, with \S n \ < 2 which is a quantitative indication 
of the rigidity of the spectrum. This is the behavior expected for a quantum chaotic system. 
The spectral fluctuations are best studied in terms of the unfolded spectrum {x n } 

x n = N(k n ) (73) 

whose mean level spacing is unity. Since ( |30[ ) provides an exact expression for N(k) the unfolding procedure is straight 
forward. 
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1. Level spacing distribution 



The distribution P(s) of the spacings s„ = x n +i ~ x n of adjacent quantal levels (or its integrated form 7(s) = 
J P(r)dr) is the most convenient and commonly used statistics. The expression for P(s) for the Poisson, GOE and 
GUE ensembles are well known, and have been compared with the distributions derived from the graph spectra. The 
numerical results for many systems show that the graphs follow the general trend observed for realistic systems. As a 
typical example, the results for the fully connected quadrangle with Neumann boundary conditions (Aj = 0) and with 
A = and A ^ are shown in Fig. 4. They are based on the leading 80,000 eigenvalues for each case. We would like 
to emphasize that the spectra were calculated for a fixed set of bond lengths, in other words, no disorder averaging 
was employed. The agreement with the exact (not the Wigner surmise) RMT curves is very good although 
systematic deviations at the level of 1% or less can be discerned (see inset of Fig. 4). These differences exceed the 
statistical error margin, and we believe that they originate from the fact that the short time dynamics on the graph 
does not follow the universal pattern, as was explained in the previous section. 

We have already noted that the Poisson limit is obtained naturally for graphs which are subject to Dirichlct 
boundary conditions. The transition between the two extremes is affected by changing the parameter A = Using 
this time a fully connected pentagon, we observe the transition in the nearest neighbor distribution as is shown in 
Fig. 5a, b for A = and A ^ respectively. 

We made similar comparisons for other well connected graphs and observed the same degree of agreement between 
the data and the results of RMT. Thus, we face an exceedingly simple class of systems which, according to the 
nearest neighbor statistics, belongs to the same spectral universality class as quantum systems which are chaotic in 
the classical limit. We shall study below other statistical measures, and show that deviations from universality appear 
as expected and observed in generic Hamiltonian systems. 



2. The form factor 



To investigate further the dynamical origins of the level fluctuations we study the two point form factor K(t; k ) 
defined by 



K(r;k ) = ^ 



E 



i 2irk n CT 



\k„-k \<A k /2 



(74) 



where we consider a spectral interval of size A*, centered about ko and involving Af = dA^ eigenvalues, r measures 
lengths in units of the Heisenberg length Ijj = C The main reason of our choice to base our analysis on the two point 
form factor is that it allows us a study of the level fluctuations in terms of PO's. Indeed, recalling that 



exp{2mkCr) = I 



fco+A fc /2 



\k n -k\<A k /2 



-Aft/2 



d(k)exp{2 / nikCT)dk 



(75) 



and expressing d(k) by its PO expansion given by ( |35| ) we can rewrite K(r;ko), after substituting the resulting 
expression into (|74|), in terms of periodic orbits. 

We shall concentrate for the time being on graphs with Neumann boundary conditions. Splitting K(r;ko) to its 
diagonal Kd(t; ko) and non-diagonal parts Knd{t; ko), we write them in terms of periodic orbits and their repetitions 
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(76) 



y ] -Apt e 



i-n(r[i p -r' n pl ) iA(rb p -r'b pl ) 
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Kb is a classical expression, because all interference effects are 



where we use A p = n p l p A p and Sj\f(x) = — ^ 

2 . . 

neglected, but for the ones which are due to exact symmetries. The sum-rule ([48|), enables us to justify a Hannay and 
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Ozorio de Almeida -like sum rule |||], namely, Kd(t) ~ (g)r p|. ((g) is the mean degeneracy of the length spectrum 
due to exact symmetries such as time reversal). For r 1 K(t) sa Kd(t). Because of the fact that the quantum 
spectrum is real and discrete, K(t) must approach 1 for r > 1. This is taken care of by Kmd- In contrast to the 
diagonal part, Knd, depends crucially on the phase correlations between the contributing terms. Actually, its Fourier 
transform tests how the l p spectrum is correlated |}9| . In Hamiltonian systems in more than one dimension, the size 
of the spectral interval is limited by the requirement that the smooth spectral density is approximately constant. 
Here d is constant, hence one can take arbitrarily large A&. This way, one can reach the domain where the function 
K(t) is composed of arbitrarily sharp spikes (Sj^(x) can become arbitrarily narrow) which resolve completely the 
length spectrum for lengths which are both smaller and larger than Lh ■ In Fig. 6a (Fig. 6b) we show the numerical 
K(t) calculated with two extreme values of Af for the case with (without) time reversal invariant symmetry. As long 
as tC is shorter than the length of the shortest period orbit, K(t) = 0, while for r > 1 it saturates and fluctuates 
around the value one. The RMT two-point form factor given as |33| 



f 2t - rln(l + 2t), < t < 1 
Kgoe ^ = \ 2 -rln(2l±i), r>l 

k gue (t) = I i' ~ > 7 1 



(77) 

is also shown in Figs. 6a and 6b for comparison. Despite the fluctuations, the low resolution curve does not deviate 
much from the prediction of RMT. The high resolution data shows a similar behavior, which can be better checked 
if one studies the integrated form factor (see inset of Figs. 6a, b) 

K(t) = - f K(t)dt. (78) 



However, by increasing the resolution, correlations between periodic orbits with different lengths are suppressed, and 
the interference mechanism which builds up K^r> cannot be due to the correlations in the spectrum of periodic orbit 
lengths, but to another source: For r > 1/2 the periodic orbits must traverse some bonds more than twice. The 
likelihood of periodic orbits which traverse the same bonds the same number of times but with different back-scatter 
indices \i p is increasing, and the interferences which build K^d are due to the sign correlations among orbits of exactly 
the same lengths (when A ^ one has to restrict the discussion to PO's with the same directed length). This result 
demonstrates one important feature of the periodic orbits correlations, namely, that periodic orbits carry not only 
metric information (lengths of trajectories) but also topological information (Maslov indices and degeneracies). The 
distribution of back scatter indices of periodic orbits is a problem that was not yet addressed by probabilistic graph 
theory. Our numerical results together with the general experience from quantum chaos allows us to conjecture that 
the spectral form factor connects RMT with the distribution of back scatter indices on PO's. 

Finally, the structure observed in the function K(r), decorating the rather smooth background can be attributed at 
low r to the short and rather scarce PO. The arrows in Figs. 6a,b indicate their location. The structures near r = 1 
reproduce a trend which was predicted on different grounds in j36fl , namely, the spikes appear at lengths C — l p (see 
arrows in Figs. 6a, b). We can explain this phenomenon with the help of Newton's identities which relate tr(S(k)) n 



to the coefficients of the characteristic polynomial, and the inversive symmetry of the latter (sec (38)-(|42|)). Simple 
algebra gives 

Y— = e* Ck+<t, ° Y ^£J_ +... (79) 

^ 2B-n ^ n K ' 

n— 1 n=l 

where the phase <po is independent of k and • • • stands for terms which involve amplitudes and phases of composite 
orbits. Substituting (|32j) and taking the Fourier transform, we find that the contributions of the terms tr5 2B_ ™ to 
the length spectrum appear at lengths C — l p where l p are lengths associated with the shorter periodic orbits with 
periods n. 

When graphs with mixed boundary conditions are investigated, (A ^ 0), the parameter A which controls the 
spectral properties depends on the mean fco parameter, and a transition from Poisson to RMT statistics is expected 
as ko increases (A decreases). This transition is illustrated in Figs. 7a,b, where the dependence of K(r;ko) on A is 
displayed. 
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3. Parametric statistics 



So far, we have shown that quantized graphs display most of the generic statistical properties encountered in the 
study of "quantum chaos". We shall discuss now yet another statistics - the parametric statistics - and show that 
the analogy carries over also for this case. Parametric statistics are defined for systems which depend on an external 
parameter (to be denoted by x), and they give a quantitative measure for the fluctuations due to level dynamics 
p4j . Among the first parametric properties studied, were the velocity distribution V(v) |Q (distribution of the first 
derivative of the levels), and the curvature distribution V(c) (distribution of the second derivative of the levels) fl46| . 
It has been shown that parametric statistics are universal for disordered or strongly chaotic systems, provided the 
change of x does not modify global symmetries. As is usually the case, non generic classical features may introduce 
deviations from the universal parametric statistics. 

The parameter which was used to study level dynamics on the graphs was the bond length of an arbitrarily chosen 
pair of bonds: Lij(x) = ^ij(O) — X an d Liiji(x) = Li'j>(0) + \, so that the total length C is kept constant. In 
this way, the mean density d is independent of x- Moreover, contrary to the usual studies of parametric statistics, 
the underlying classical dynamics of the graph are unaffected by the change of x ( see also @ ) ■ Modulating the two 
other parameters of our graphs i.e. the "magnetic potential" A and the scattering potential A at each vertex, we are 
able to study the parametric statistics in the transition regions where either time reversal symmetry or integrability 
are broken. 

To .evea, the universality in „(„) and V(c), one «e a the variance a, = to tescaie the verities . and 

the curvatures c 



k- 



dk n d 2 k„ 
dx „ _ dx 2 
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(80) 



where /3 = 1, 2 correspond to graphs with or without time reversal symmetry, respectively. 

The numerical analysis reported below was conducted on a fully connected hexagon. For intermediate boundary 
conditions the wavenumber range, over which we performed our statistical analysis was kept small in order to keep the 
control parameter A to be essentially constant. Then the statistics was generated over realizations of the lengths of the 
bonds. By employing a finite-difference method, we were able to compute the level velocities and the curvatures for 
many different values of the parameter x- The total number of eigenvalues used to construct the histograms exceeded 
186000 in each statistics. 

We first analyze the velocity distributions for graphs. For level dynamics within the GOE or the GUE the distri- 
bution of level velocities V{v) is proved to be Gaussian [H. Some of our results are shown in Fig. 8a, b for A = 
and A ^ 0, respectively (in both cases, Neumann boundary conditions where imposed). The calculated velocity 
distributions are well approximated by a Gaussian of the same mean value and standard deviation. 

Fyodorov ^8| derived an analytical formula for V(v) which applies for disordered systems in the strongly localized 
limit and for which time reversal symmetry is violated 

■n irv coth (ttv I v6) — 
~ 6 smnVw/v^) ' 

Numerical simulations have shown that this formula is also applicable in the domain where the Poisson-GOE transition 
takes place ^7j. Our numerical results, presented in Fig. 9 show that (|8l]) reproduce the data in the range of large 
A values, thus confirming the suggestion first made in p7[ . 

RMT predicts explicit expressions for the curvature distributions V(c) |Q 

^>- % + c2 )(. +2 )/ 2 ( 82 ) 

where Np equals to 0.5 and 2/tt for GOE ((3 — 1) and GUE {(3 = 2) respectively. Our numerical calculations for the 
cases with Neumann boundary conditions, and A = 0, A ^ are shown in Figs. 10a,b respectively. The agreement 
with the theoretical expectation ( |S2| ) is excellent. 

In Ref. jl7|, it was suggested that d82| ) can be generalized for the intermediate statistics interpolating between 
Poisson and GOE or GUE. The normalization constant Np has to be defined as 



^r( 2 



N ? - — ^TTT- ( 83 ) 
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and (3 takes real values within the interval (0, 2] (see |49]). We checked this conjecture for our system and we found 
that it describes in a satisfactory way the intermediate statistics. For this, we fitted the tails of the distribution of 
the second derivative of the levels (unsealed curvatures) to an algebraic decay V(G) ~ C~^ +2 K Then the value of (5 
found from the fit was used in ( |80| ) to rescale the curvatures. Our results for the case of partially broken time reversal 
symmetry with A = 4 are shown in Fig. 11a. Similarly, in Fig. lib we report our findings for the transition between 
Poisson and GOE. 



VI. GRAPHS WITH NON UNIFORM CONNECTIVITY 



So far, we have studied properties of well connected graphs, and have shown that when the appropriate limit is taken, 
many statistical properties of the spectrum reproduce the expectations of RMT. We shall dedicate the next section 
to demonstrate cases for which the connectivity of the graph induces non uniform dynamics which has substantial 
effects on the corresponding quantum spectra and their statistics. 



A. The Hydra 

As a first example of a family of graphs which is not uniformly connected, we studied the "Hydras" or "star" graphs. 
They are graphs which consist of v$ bonds, all of which emanate from a single common vertex labeled with the index 
i = 0. The vertex at i = 0, will be referred to in the sequel as the Hydra's head. The total number of vertices for 
such a graph is V = Vo + 1, and the vertices at the end of the bonds will be labeled by i = 1 • • • V. We shall assume 
Neumann boundary conditions on these vertices. The Hydra is a bipartite graph, a property which implies e.g., that 
there exist no periodic orbits of odd period! This is responsible for most of the non generic properties of the classical 
and the quantum properties of Hydras. 

We start with the Sb~ matrix statistics which will allow us a better physical understanding of our system. To this 



end we first calculate the matrix T defined in (17). One can easily show that 



T =(° (0) o)- (84) 

The matrix is the vq x vq scattering matrix at the Hydra's head as defined in (pd|). / denotes the vq x vq identity 
matrix. It represents the trivial back scattering at the vertices i = 1 • • • V. The Sb (k; A) matrix and its square 
S% (k\ j4)can be written as 

S B {k;A) = y d (- )ai0 ) Q j 5 S B (k;A)=^ Q d(-) CT (o) d (+) ) ( 85 ) 

Where ^ — 5i.je l ( k±A ^ Lj are the diagonal matrices which carry the metric information. It follows from ( p5| ) that 
trS'^ n+1 = 0, and the even traces satisfy trS^" = 2trS'^, where Sh = d^a^°'d^ + '. Sh is a v x vq scattering matrix 
in the space of the Hydra bonds. It incorporates the reflections from the vertices i = 1 • • • V. Sh is independent of A 
since all the PO's on the Hydra are self tracing, and from now on we shall study its properties, since it is free from 
the trivial effects which originate from the bipartite nature of the Hydra. 

If we consider the limit vq 3> 1, and for n < vq, we can use the diagonal approximation to calculate the form factor 
of the Sh matrix eigenphase spectrum. In the present context, the classical evolution operator which corresponds to 
the quantum Sh matrix is 

(^kj = ^C 1 + COSUJ o) + cosw o)) Ko ■ (86) 

The spectrum of Uh consists of the values 1 and 1 — 2 1 + c ° s "f which is Vq — 1 times degenerate. Therefore, 

u n = txU H = l + (v -l)(l- ^-(1 + coscjo)^) • (87) 

For large valencies vq the t?S H are dominated by the n repetitions of fixed points. Thus, we can write for the trace 
oiS H 
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trS H = "£e 
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where the second sum contains contributions from other periodic orbits. By performing furthermore an average over 
realizations of the lengths of the bonds of Hydras we get 
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(89) 



In the second term of the r.h.s of ( |89| ) we have subtracted from the classical return probabilities u n (see (J8J), the 
contributions from the fixed points which already had been taken into account. ( |89| ) can be written in a better form 
i.e. 
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(90) 
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This is different from the generic expression in two important ways: The linear term which dominates the larger n 
domain is proportional to n and not to In because in the Hydra all the periodic orbits are self tracing. The other two 
terms which dominates the domain of smaller n, are due to the n repetitions of fixed points and to small (degenerate) 
eigenvalues of Ur- Their n dependence cannot be scaled with vq, and they represents a typical transient effect. The 
spectrum of Sh is not degenerate, and therefore for n > vo the form factor reaches its asymptotic value 1. This 
transition is due to the interference of contributions of periodic orbits with the same lengths, but even for the simple 
Hydra graphs we do not have a theory which explains this phenomenon. 

The above approximation applies quite well also for the Dirichlet boundary conditions. Indeed, for loq = tt we have 
from j9C| ) that ^ (\teS H \ 2 ^ = 1. Moreover, the transition to the Poisson limit is again described quite well by (|9(] ), 
for r = ^- < 1. This is shown in Fig. 12 for a Hydra with v = 50 and various values of the parameter A . The non 
generic features expected from the periodic orbit theory are well reproduced by the numerical data. 

The wavenumber spectrum follows the same trends as the eigenphase spectrum discussed above. The secular 
equation (^|) for Hydras takes a particularly simple form, 
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Here, we investigate for simplicity the case of Hydras with a zero scattering potential Ao = 0. In Fig. 13 we present 
the level spacing distribution P(s) for various valencies. In all cases, the spectrum shows level repulsion, which is 
described rather well by the Wigner surmise for large valencies. A closer look shows deviations which do not decrease 
when the valency increases. In particular, in the limit s — > 0, the 7(s) does not approach the expected power law, and 
this can be seen in the inset of Fig. 13 where the results for Hydras with vq = 7 — 18 are compared with the GOE 
expression. In light of the discussion of the Sh eigenvalues statistics, the deviations observed in the nearest neighbor 
spacings distribution are not surprising. 

In a previous publication ]23| we investigated the two point form factor K(t) for the Hydras and we had reported 
(for the case Vq = 5) quite good agreement with RMT predictions. However, upon increasing Vo, K(t) shows the same 
behavior as its counterpart for the eigenphase spectrum of Sh- In Fig. 14 we present the two-point form factor for 
Hydras with ^o = 5, 15. The deviations from RMT predictions, as we are increasing the valency, follow the pattern 
observed for the Sh spectral statistics (see (|90|)). 

As expected, the Hydra spectrum displays non generic spectral statistics, which reflect the special connectivity 
and hence the classical dynamics of this graph. Periodic orbits theory (in particular, its diagonal approximation) 
reproduces the statistics which relate to the limit of short evolution times. 



B. Rings With Variable Connectivity 

In this section we shall study a family of regular graphs which are defined by the connectivity matrix 
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_ J 1 for i ^ j and \i - j\ < b , . 

lJ ~ \ for i = j or |i - j\ > b {y > 

so that 2b is the valency. When the graph is drawn in the plane, with the V vertices placed on a circle with 
9i > 8j for i > j, and V ^> b, then the resulting shape is a ring, and hence the name. 

We shall study the spectral statistics as a function of the valency, and will attempt to answer the question, at what 
value of b can we consider the graph as sufficiently well connected so that its spectral properties can be reproduced 
by RMT. We shall assume Neumann boundary conditions throughout this section. 

As long as 1 < b <C V, the classical evolution operator describes a random walker on a ring, where the hopping 
step size can take any value between 1 and b. This results in a diffusive evolution on the ring, with D w f tt 2 6 2 . The 
case b = 1 is trivial because the Neumann boundary conditions in this case does not allow reflections at the vertices, 
and hence, the particle goes around the ring "ballistically" . To eliminate this effect we added loops at each vertex. 
This does not alter the diffusive dynamics, but allows us to include rings with nearest neighbor hopping. 

The quantum dynamics is strongly affected by the fact that the lengths of the bonds are rationally independent. 
In the limit b <C V this causes the eigenstates to be localized, and the spectral correlations will bear the marks of the 
degree of localization. Our purpose is to study the transition in the spectral statistics of rings, as b is increased. The 
addition of loops to the vertices, introduces also a subspace of eigenfunctions for which the vertex values <f>i = 0, and 
the wavefunction itself vanishes on all the bonds and all the loops but one. There it takes the form sin jttjX, where 

Lf is the length of the loop attached at the i vertex. This subspace of eigenstates and the corresponding spectrum 
will be excluded from the following statistical study. 

Starting with rings with 6=1, we can calculate their properties by introducing the quantum transfer operator 

' Ti -[ P (k,A^) l-MMV AO eW^J (93) 

where the superscript (I) distinguishes the loops parameters. The bond (i — is referred to as the i bond. The 
parameter p is given by 

. A (l) +k T U) . A ll) -k T (l) 

,f;v sin— Hj — L,- sin — «-= — L) ' , 

p(k, Af ) = 2i 2 - ! ^ l - . (94) 

sin kL I ' 

By multiplying the transfer matrices for arbitrarily long segments, we calculated the Lyapunov exponent which 
yields the localization length Z^. Since for all k values we found it to be of order 1, we expect that the eigenfunctions 
for rings with V > 10 will be well localized. A direct inspection of the eigenstates confirms this expectation (see 
Fig. 15). 

The effect of localization on the spectral statistics can be understood by the following argument. As long as the 
total length C of the ring is sufficiently larger than the Anderson localization length Zoo, one can approximate the 
spectrum qualitatively as a union of Cllca uncorrelated spectra. One expects that the form factor for two rings with 
lengths C and C will be related by |5(j] 

K c {t)^Kc{tC/L'). (95) 

In Fig. 16 we plot the integrated form factor K(t) (|t^) for b = 1 graphs with number of vertices V — 11 and V = 22. 
Scaling the r axis by 2, the two form factors coincide as expected. The relation (|9^) provides the correct sense in 
which one should interpret the statement that the spectrum of a localized system tends to the Poissonian limit as the 
size increases. The correlations in the spectrum remain, but they are on a scale which is C/lco larger than the mean 
spacing. 

As the connectivity is increased, and keeping the length constant, one expects to see a transition in the spectral 
statistics. The first statistical quantity that we investigate is the integrated level spacing distribution I{s). Our 
results for the case V = 22 are shown in Fig. 17. As the connectivity range increases, I(s) makes a transition from 
the Poisson distribution towards the GOE. Recently, it was suggested |3l| that there exists a universal intermediate 
statistics, which applies for systems which undergo a transition from Poisson to RMT statistics. At the intermediate 
(critical) point, 

P cr (s) = 4sexp(-2s). (96) 
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Fig. 17 shows an excellent agreement between our numerical results for the ring with b = 2 and (|96|). Similar degree of 
agreement with the critical statistics appears also for other statistical measures like the number variance. The number 
variance S 2 (L) probes the spectrum over all correlation lengths L and describes the fluctuations of the number n(L) 
of levels contained in a randomly chosen interval of length L. Is defined as 

T?{L)^{{n{L)-Lf) (97) 

where the angular brackets () denote a local averaging over sufficiently many levels. S 2 (L) is related with K(t) by 
the integral transform 

^Lj^frfr^M). (98) 



and the reason that we choose here to concentrate on S 2 (L) is that it behaves quite nice, with respect to the highly 
fluctuative form factor. 

The general expectation for generic systems p3| ], is that £ 2 (L) should comply with the predictions of RMT for small 
values of L (universal regime) and saturate to a non- universal value for large L's due to the semiclassical contributions 
of short periodic orbits. Our results are presented in Fig. 18 together with the critical number variance which reads 

13 

Zl(L) = l + 1 - eXP 8 { - 4L) . (99) 

For reference we also draw the number variance for a Poissonian spectrum given by ^(L) = L, and the RMT 
expectation 



£ 2 (£) GO b= 4{ ln ( 27ri ) +7+ 1 + Jsi 2 (7ri) - J Si(7rL) - cos(2ttL) - (100) 
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where 7 is the Euler constant. 

As one increases the connectivity the graph statistics approach the RMT expressions. On the basis of the estimate 
of the diffusion constant on the ring, we expect that the RMT limit is reached for b ~ V. 



VII. CONCLUSIONS 



The graph is a one dimensional system. Yet, it is not simply connected, and this is why it can display chaotic 
classical dynamics, in the sense explained above. The fact that the graph is one dimensional can be seen in various 
classical and quantal attributes. On the quantum level, the smooth spectral counting function is proportional in 
general to E d t 2 and indeed, in our case d = 1. The stability factors in the trace formula correspond to a system with 
a single expanding direction, and no contracting direction. In this respect the graph trace formula is similar to the 
Riemann-Weyl trace formula. (The counting function for the Riemann zeros, however, has a logarithmic correction 
to the strictly one dimensional Weyl term). 

In spite of the simplicity of the graph Schroedinger operator, we have shown that the spectrum displays many 
features which appear in the study of the quantum analogues of classically chaotic systems. The limit h — > 00 can be 
replaced by the limit C — > 00, which ensures that all the features which are due to the short, and non generic periodic 
orbits are of lesser relative importance. 

For a well connected graph, the length spectra of periodic orbits of periods larger than V show increasing degree 
of degeneracy as the period increases. These degeneracies give rise to the correct behavior of the form factor at large 
values of r. In other words, the fact that the quantum spectrum is real and discrete, is expressed in the trace formula 
through the degeneracies of periodic orbit lengths, and the distribution of the corresponding back scattering indices. 
This observation points at a possible new interpretation of the Wigner Dyson theory in terms of probabilistic graph 
theory. 

The results reported here encourage us to believe that quantum graphs may serve as a convenient tool for the study 
of quantum chaotic or disordered systems. We do not know yet if it is possible to map any given Hamiltonian system 
onto a corresponding graph. The evidence is mounting, however, in favor of a very intimate link, and the search for 
this connection is one of our immediate goals. 



24 



VIII. ACKNOWLEDGMENTS 



This research was supported by the Minerva Center for Physics of Nonlinear Systems, and by a grant from the 
Israel Science Foundation. (T. K.) thanks the Feinberg School at the WIS for a Post Doctoral fellowship. We thank 
D. Cohen, D. Miller, H. Primack, Z. Rudnick and H. Schanz for valuable suggestions and comments, and Y. Colin 
De Verdiere for bringing J. P. Roth's paper jl7| to our attention. In particular we thank U. Gavish for reading the 
manuscript critically and making valuable comments. 



25 



L. Pauling, J. Chem. Phys. 4, 673 (1936) 

H. Kuhn, Helv. Chim. Acta, 31, 1441 (1948); H. Kuhn ibid, 32, 2247 (1949). 
J. R. Piatt, J. Chem. Phys. 17, 484 (1949). 

K. Rudenberg, J. Chem. Phys. 21, 1509 (1952); K. Rudenberg and C. W. Scherr, J. Chem. Phys. 21, 1565 (1953). 
C. A. Coulson, Proc. Phys. Soc. Lond. 67, 608 (1954). 
E. W. Montroll, J. M. P. 11, 635 (1970). 

Ph. Choquard, Lecture notes Phonon end Electrons II, section 3.1 (1971). 
M. J. Richardson, N. L. Balazs, Annals of Physics 73, 308 (1972). 
S. Alexander, Phys. Rev. B 27, 1541 (1985). 

C. Flesia, R. Johnston and H. Kunz, Europhys. Lett. 3 , 497 (1987). 

R. Mitra and S. W. Lee, Analytical techniques in the theory of guided waves, Macmillan, New York. 
P. W. Anderson, Phys. Rev. B 23, 4828 (1981); B. Shapiro, Phys. Rev. Lett. 48 823 (1982). 

J. T. Chalker and P. D. Coddington, J. Phys. C 21, 2665 (1988); Rochus Klesse and Marcus Metzler, Phys. Rev. Lett. 79, 
721 (1997); R. Klesse, Ph. D. Thesis, Universitat zu Koln, AWOS-Verlag, Erfurt, 1996. 

Y. Avishai, J. M. Luck, Phys. Rev. B 45, 1074 (1992); T. Nakayama, K. Yakubo, R. L. Orbach, Rev. Mod. Phys. 66, 381 
(1994). 

Y. Imry, Introduction to Mesoscopic Systems, Oxford, New York (1996); D. Kowal, U. Sivan, O. Entin-Wohlman, Y. Imry, 
Phys. Rev. B 42, 9009 (1990). 

J. E. Avron, in: Proc. 1994 Les Houches Summer School on Mesoscopic Quantum Physics, E. Akkermans et.al., eds. 
(North-Holland) 741-791. 

Jean-Pierre Roth, in: Lectures Notes in Mathematics: Theorie du Potentiel, A. Dold and B. Eckmann, eds. (Springer- 
Verlag) 521-539. 

P. Exner, P. Seba, Rep. Math. Phys. 28, 7 (1989); P. Exner, J. Phys. A: Math. Gen. 29, 87 (1996). 
P. Exner, R. Gawlista, Phys. Rev. B 53, 7275 (1996); P. Exner, Phys. Rev. Lett. 74, 3503 (1995) 
J. E. Avron, P. Exner, Y. Last, Phys. Rev. Lett. 72, 896 (1994). 

R. Carlson, "Inverse Eigenvalue Problems on Directed Graphs", Transactions of the American Mathematical Society, to 
appear; R. Carlson, Electronic Journal of Differential Equations 23, 1 (1997); R. Carlson, Electronic Journal of Differential 
Equations 6, 1 (1998). 

V. Kostrykin and R. Schrader, Kirchhoff's Rule for Quantum Wire, preprint. 
T. Kottos and U. Smilansky, Phys. Rev. Lett. 79, 4794 (1997). 

M. G. Gutzwiller, in Chaos in Classical and Quantum Mechanics, Vol. 1 of Interdisciplinary Applied Mathematics, edited 
by F. John (Springer- Verlag, New York, 1990). 
J. P. Keating, Nonlinearity 4, 309 (1991). 

U. Smilansky, in Proc. 1994 Les Houches Summer School on Mesoscopic Quantum Physics, E. Akkermans et.al., eds. 
(North-Holland) 373-433. 

J. Gratus, C. J. Lambert, S. J. Robinson and R. W. Tucker, J. Phys. A: Math. Gen. 27, 6881 (1994). 
T. Kottos and U. Smilansky, in preparation. 
A. Selberg, J. Ind. Math. Soc. 20, 47 (1956). 

E. Doron and U. Smilansky, Nonlinearity 5, 1055 (1992). 

J. Keating and M.V. Berry, Proc. Royal Soc. Lond. A 437, 151 (1992). 

U. Smilansky, in Proc. 1989 Les Houches Summer School on Chaos and Quantum Physics, M.-J. Giannoni et.al., eds. 
(North-Holland) 371-441. 

M. L. Mehta, Random Matrices and the Statistical Theory of Energy Levels, Academic, New York (1990); T. A. Brody, J. 
flores, J. B. French, P. A. Mello, A. Pandey, S. S. M. Wong, Rev. Mod. Phys. 53, 385 (1981). 

R. Blumel and U. Smilansky, Phys. Rev. Lett. 60, 472 (1988); R. Blumel and U. Smilansky, Phys. Rev. Lett. 64, 241 
(1990); 

F. Haake, M. Kus, H. -J. Sommers, H. Schomerus, and K. Zyckowski, J. Phys. A 29, 3641 (1996). 

E. B. Bogomolny and J. P. Keating, Phys. Rev. Lett. 77, 1472 (1996); O. Agam et.al., ibid 75, 4389 (1995). 
U. Smilansky, Physica D 109, 1767 (1997). 

U. Smilansky, Semiclassical Quantization of Maps and Spectral Ccorrelations, Proc. 1997 of the Nato Advanced Study 
Institute "Super symmetry and Trace Formulae ", I. Lerner, eds. (Cambridge) (in press). 

N. Argaman, F. M. Dittes, E. Doron, S. P. Keating, A. Y. Kitaev, M. Sieber and U. Smilansky, Phys. Rev. Lett. 71, 4326 
(1993) ; D. Cohen, H. Primack and U. Smilansky, Annals of Physics 264, 108-170, (1998). 
S. Kettemann, D. Klakow and U. Smilansky, J. Phys. A: Math. Gen. 30, 3643 (1997). 
Barbara Dietz, F. Haake, Z. Phys. B 80, 153 (1990). 

J. H. Hannay and A. M. Ozorio De Almeida, J. Phys. A: Math. Gen. 17, 3429 (1984). 



2G 



[43] M. V. Berry, Proc. Royal Soc. Lond. A 400 229 (1985). 

[44] J. Goldberg, U. Smilansky, M. V. Berry, W. Schweizer, G. Wunner and G. Zeller, Nonlinearity 4, 1 (1991). 
[45] B. D. Simons, B. L. Al'tshuler, Phys. Rev. B 48, 5422 (1993); 

[46] P. Gaspard, S. A. Rice, H. J. Mikeska, K. Nakamura, Phys. Rev. A 42, 4015 (1990); J. Zakrzewski, D. Delande, Phys. 

Rev. E 47 , 1650 (1993); F. von Oppen, Phys. Rev. E 51, 2647 (1995); F. von Oppen Phys. Rev. Lett. 73, 798 (1994). 
[47] P. Kunstman, K. Zyczkowski, J. Zakrzewski, Phys. Pev. E 55, 2446 (1997). 

[48] Y. V. Fyodorov, Phys. Rev. Lett. 73, 2688 (1994); Y. V. Fyodorov, H.-J. Sommers, Phys. Rev. E 51, R2719 (1995); Y. V. 

Fyodorov, H.-J. Sommers, Z. Phys. B 99, 123 (1995). 
[49] F. M. Izrailev, Physics Reports 129, 299 (1990). 
[50] Pandey A., Ann. Phys. (NY) 134 110 (1981). 
[51] E. B. Bogomonly, in preparation. 



27 




1 

"T 












i.i. 



0.5 

— i 



.O 

-I— 



.5 

— i — 




O.O 0.5 1.0 1.5 



1 .5 



n/2B 



2.0 



Fig. 1: The form factor for the eigen-phase spectrum of 5s for a pentagon graph. Bold solid lines are the expecta- 
tions for the COE and CUE expressions. The data are averaged over odd-even powers of the form factor as explained 
in the text. In the insets we present the form factor for the case A = 0.833 (o) without averaging. 

(a) A = 0, A = (□), A = 1.25 (*), A = 2.5 (o)and A = 3500 (•). 

(b) A ^ with the same boundary conditions as in (a) 
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Fig. 2: The Fourier components (|a;| 2 ) of Cq for a pentagon graph. (*) corresponds to semi-classical calculations 
while (o) to the exact quantum mechanical calculations. 

(a) A = 0, A = 25000, A = 25, A = 0.476 and A = 0; 

(b) A ^ 0, with the same boundary conditions. 

In both figures the bold solid line correspond to Poisson while the bold dashed line is the expectation of the corre- 
sponding RMT. 
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Fig. 3: S n = N(k) — N(k) vs. the wavenumber label n for the graph shown in the inset. 
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Fig. 4: Integrated nearest neighbor distribution I(s) for a fully connected quadrangle with A = (Neumann bound- 
ary conditions). The results are based on the lowest 80, 000 levels of a single realization of the bonds. AI indicates 
the deviation from the exact GOE/GUE results; 
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Fig. 5: The integrated nearest neighbor distribution I(s) for a fully connected pentagon and various values of the 
parameter A. The statistics was generated over a large number of realizations of the bonds of the graph; 
(a) A = and various values of the parameter A. (b) A ^ with the same boundary conditions as in (b). 
In both figures the thick solid line is the expectation of the corresponding RMT while the thick dashed line correspond 
to Poisson . 



30 



1 O EE ■ 1 ■ 1 ■ 1 ■ 1 ■ r 



K( 



K(x 




Fig. 6: Two point form factor (w 100, 000 levels). The arrows indicates the location of the short periodic orbits and 
their reciprocal lengths with respect to the Heisenberg length. In the insets we show the corresponding integrated 
form factor K(t) (thin line). 

(a) Fully connected pentagon with ^4 = 0; 

(b) Fully connected quadrangle and A ^ 0. 

In both figures the bold lines are the expectation of the corresponding RMT. 
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Fig. 7: The integrated two point form factor K(t) for a fully connected pentagon. The thick dashed lines correspond 
to the RMT expectations. 

(a) A — and various values of A. 

(b) A 7^ and various values of A. 
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Fig. 8: Velocity distribution P(v) for a fully connected hexagon. The thick dashed line correspond to a Gaussian 
with the same mean and standard deviation, with the numerical data; 

(a) A — 0, Aj = (Neumann boundary conditions). 

(b) A + 0, A, = 0. 
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Fig. 9: Velocity distribution P(v) for a fully connected hexagon. A = 0, and Aj = A = 4.97. The thick solid line 
correspond to ( |8l| ) which describes the Poisson-GOE transition. The thin solid line correspond to a Gaussian with 
the same mean and standard deviation, with the numerical data. 
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Fig. 10: Curvature distribution P(c) for the fully connected hexagon. The thick dashed line correspond to (|8^ ) 
with P defined from the symmetry of the system. 

(a) A — and A, = (Neumann boundary conditions). The parameter j3 = 1. 

(b) A = 20 and A* = 0. The parameter /3 = 2; 
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Fig. 11: Curvature distribution P(c) for the fully connected hexagon. The thick dashed line correspond to ( p^ ) 
with (3 defined from the symmetry of the system (fitting parameter). 

(a) A — 2 and A, = 0. The fitting parameter is f3 = 1.4 ± 0.1; 

(b) A = and Aj = A = 4.97. The statistics was generated as explained in the text. The fitting parameter is 
[3 = 0.28 ±0.01. 




Fig. 12: The phase-shift form factor for a Hydra with u = 50 and A = 2000 (o), A = 2 (□) and A = (O). The 
thin solid lines correspond to the theoretical expectation (|90|), while the thick dashed line to the RMT prediction. 
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Fig. 13: Nearest level spacing distribution P(s) for various valencies vo of the Hydras and Ao = 0. In the inset we 
report the integrated level spacing J(s) vs. s, for small spacings (the dashed lines are guides for the eye). Deviations 
from RMT (thick solid line) are clearly observed. 




Fig. 14: The two-point form factor K(t) for Hydras with vo = 5 (thin solid line) and vq = 15 (thin dashed line). 
The thick solid line correspond to the RMT expectations. In the inset we show the integrated form factor K(t). 
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Fig. 15: A typical exponentially localized eigenstate of a graph with V = 44 and loops at each vertex. 
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Fig. 16: The integrated form factor K(t) for graphs with V = 11 (A) and V = 22 (□) and loops at each vertex. 
The thick solid line correspond to the graph with V — 22 after rescaling r by a factor of two. 
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Fig. 17: Integrated nearest neighbors distribution J(s) for a graph with loops at each vertex. The number of vertices 
is V — 22. We consider various connectivities; 6=1 (□), 6 = 2 (o), 6 = 3 (*). The critical distribution ( p6| ) is also 
shown (thick dashed line) together with the GOE (thin solid line) and Poisson (thin dashed line) expressions. 
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Fig. 18: The number variance S 2 (L) for the graph of Fig. 17 and various connectivities; 6=1 (□), 6 = 2 (o), 6 = 3 
(*), 6 = 5 (A). The critical distribution ( p^ ) is shown by the thick dashed line, together with the GOE (thin solid 
line) and Poisson (thin dashed line) results. 



37 



